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ABSTRACT 


INTERACTION  OF  ELECTROMAGNETIC 
SIGNALS  WITH  PARTICULATE 
CLOUDS 

by 

JOSE  MARIO  PAUDA,  B.S. 

SUPERVISING  PROFESSOR:  THOMAS  A.  GRIFFY 

A  particulate  cloud  affects  the  ability  of  an  electronic  detector  to 
receive  an  electromagnetic  signal  in  two  ways:  by  scattering  light  from  the  sun 
into  the  detector,  thereby  masking  the  signal,  and  by  attenuating  the  signal 
itself.  These  effects  are  well  studied  in  the  Mie  theory,  which  is  summarized. 
The  effect  of  the  particle  distribution  in  the  cloud  and  the  shape  of  the  cloud 
on  the  scattering  and  absorption  problems  is  then  analyzed.  The  results  of 
this  analysis  and  of  the  Mie  theory  are  incorporated  into  a  computer  program 
which  is  included  in  the  appendix.  The  graphs  generated  with  the  program  can 
be  used  (in  conjunction  with  information  about  the  sunlight  intensity  and  the 
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Chapter  1 


Introduction 

1.1  Scattering:  A  New  Computer  Model 
for  an  Old  Problem 

Most  of  the  information  we  receive  from  a  variety  of  sources  is  in  the 
form  of  electromagnetic  signals.  These  electromagnetic  signals  interact  with 
the  particles  suspended  in  the  medium  through  which  they  travel.  Examples 
include  the  absorption  of  sunlight  by  dust  particles  suspended  in  the  atmo¬ 
sphere,  the  absorption  of  radar  signals  by  water  dropplets,  and  the  absorption 
of  infrared  signals  by  particles  produced  by  rocket  motors.  If  we  assume  that 
an  electronic  detector  is  to  receive  and  process  a  signal  that  has  to  interact  with 
a  cloud,  be  it  a  radar  return  passing  through  a  rain  cloud  or  an  infrared  signa¬ 
ture  of  a  “hot”  object  passing  through  a  smoke  cloud,  we  must  investigate  the 
effect  of  the  particulate  cloud  on  the  signal  itself.  This  effect  is  two-fold:  the 
signal  may  be  partialy  absorbed  by  the  particles  to  a  point  below  the  thresh¬ 
old  of  the  detector,  or  the  signal  may  be  masked  by  the  light  from  a  brighter 
source  scattered  by  the  particles  into  the  detector.  This  is  the  major  aim  of 
this  project:  to  create  a  computer  model  to  investigate  the  absorption  (also 
called  the  extinction)  of  the  signal  as  it  travels  through  a  cloud  of  particles  and 
the  masking  of  the  signal  by  light  scattered  by  a  cloud  of  particles.  Specifically, 
we  consider  the  effect  a  conical  cloud  has  on  a  detector  at  the  vertex  of  the 
cloud.  This  investigation  is  by  no  means  new.  However,  other  authors  have 
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limited  their  investigations  to  long  wavelengths  or  to  preselected  absorption 
values.  We  have  tried  to  fill  this  gap  by  providing  a  general  computer  code, 
with  none  of  the  limitations  listed  above.  Thus,  by  implementing  the  code 
included  in  the  appendix,  the  reader  will  be  able  to  generate  theoretical  pre¬ 
dictions  for  a  wide  range  of  particles,  wavelengths  and  refraction  coefficients. 
In  the  remainder  of  this  chapter,  a  brief  explanation  of  how  these  parameters 
affect  the  absorption  and  scattering  of  the  signal  is  presented.  In  particular, 
the  results  of  the  Mie  Theory  are  summarized  and  expressions  for  the  electric 
and  magnetic  fields  in  terms  of  known  parameters  are  presented.  In  chapter  2 
the  effect  a  nonconstant  particle  size  distribution  is  investigated,  as  well  as  the 
effect  of  the  shape  of  the  cloud  itself.  In  chapter  3  we  apply  the  information 
from  chapters  I  and  2  and  present  the  graphs  of  the  scattered  light  and  ab¬ 
sorbed  signal.  These  graphs  were  obtained  by  plotting  the  points  calculated 
by  the  computer  program  included  in  the  appendix. 

1.2  Exact  Solution  to  the  Problem: 

The  Mie  Theory 

Any  electromagnetic  signal  may  be  treated  as  a  plane  wave  if  the 
distance  from  the  source  to  the  observation  point  is  much  greater  than  the 
wavelength.  More  explicitly,  if  the  observation  point  is  at  a  distance  d  from 
the  source  and  the  signal  has  wavelenght  A,  then  with  k  =  the  signal  may 
be  treated  as  a  plane  wave  if 

kd>  1.  (1.1) 

In  particular,  consider  a  spherical  particle  in  a  vacuum  at  a  distance  d  from 
the  sun,  such  that  equation  1.1  holds.  Furthermore,  assume  that  the  particle  is 
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Figure  1.1:  Relation  of  the  incoming  plane  wave  to  the  scattering  angle  0 
[van  de  Hulst  81,  page  12]. 


absorbing  in  the  ..lfrared  region  of  the  spectrum.  Then  the  scattering  problem 
is  as  depicted  in  figure  1.1:  a  train  of  plain  waves  comes  in  from  “infinity” 
and  interacts  with  the  particle.  Part  of  the  signal  is  absorbed  and  part  is 
scattered  at  an  angle  6.  If  the  radius  of  the  particle  is  much  smaller  than  the 
wavelength,  r<l,  the  problem  reduces  to  the  well-known  Rayleigh  scattering, 
treated  in  many  textbooks.  Jackson,  for  instance,  has  several  sections  dedicated 
to  this  very  problem  [Jackson  75,  sections  9.6  and  9.7].  However,  to  solve  the 
scattering  problem  exactly,  with  no  assumptions  as  to  the  ratio  of  the  particle’s 
size  to  the  wavelength  of  the  signal,  we  must  solve  Maxwell’s  equations  subject 
to  the  appropriate  boundary  conditions.  This  problem  was  first  solved  by  Mie 
in  1908  [Mie  08].  In  his  classic  book  van  de  Hulst  covers  the  work  of  Mie 
thoroughly  and  concisely  [van  de  Hulst  81,  chapter  9].  Indeed,  his  treatment 
of  the  scattering  of  light  is  among  the  most  lucid  and  complete.  In  this  thesis, 
we  do  not  recreate  Mie’s  original  work,  nor  emulate  van  de  Hulst’s.  Instead,  we 
simply  summarize  the  process  they  used  to  arrive  at  the  solution.  Their  result 


4 


is  then  used  to  write  the  computer  program  which  can  be  used  to  calculate  the 
effects  of  a  cloud  of  particles  on  the  signal. 

Consider,  then,  a  plane  wave  of  random  polarization,  and  a  spherical 
particle  with  which  the  wave  interacts,  as  in  figure  1.1.  Let  the  electric  field  be 
represented  by  E  and  the  magnetic  field  be  represented  by  B.  Furthermore, 
let  the  wave  have  wavelength  A  and  circular  frecuency  ui  =  ck.  Then,  in  any 
medium,  this  wave  is  represented  by  Maxwell’s  equations  as1 

V.£  =  0  V  x  £  +  7^=  0 
V-B  =  0  V  X  =  0 

where  c  is  the  speed  of  light  in  vacuum,  /i  is  the  magnetic  permeability  and 
e  is  the  electric  permittivity  of  the  medium  (in  vacuum,  these  parameters  are 
unity).  By  combining  the  two  curl  equations,  we  can  show  that  each  cartesian 
component  of  E  and  B  satisfies  the  wave  equation 

V’u-i|H  =  0  ,1.2) 

where 

»  =  -f-  (1-3) 

y/Ht 

is  the  velocity  of  the  wave  in  the  medium.  For  a  wave  traveling  in  the  z 
direction,  the  wave  equation  1.2  has  the  solution 


u(z,  t)  =  e 


_  — «fcz+*u/t 


lIn  this  introductory  chapter,  we  follow  the  approach  of  van  de  Hulst,  chapter  9.  In 
particular,  we  use  his  notation  for  a  plane  wave  traveling  in  the  z  direction  as  shown  in 
equation  1.4.  Most  modern  references  assign  a  negative  sign  to  the  time  component  of  the 
wave  and  a  positive  sign  to  the  spacial  component,  which  is  opposite  to  the  convention 
adopted  here.  In  addition,  we  have  used  the  cgs  unit  system. 
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This  solution  may  be  generalized  to  include  the  vector  quality  of  the  E  and  B 
fields.  Following  the  convention  that  the  physical  electric  and  magnetic  fields 
are  obtained  by  taking  the  real  parts  of  complex  quantities,  the  plane  wave 
fields  may  be  written  in  the  form 


E(z,t)  =  e\Ee~ik2+iwt 

(1.5) 

B(z,  t)  =  e2Be~ikz+iwt 

(1.6) 

Equations  1.5  and  1.6  describe  a  plane,  transverse  wave  in  vacuum,  with  the 
electric  field,  of  magnitude  E,  polarized  in  the  ej  direction  and  the  magnetic 
field,  of  magnitude  B,  pointing  in  the  e2  direction.  The  waves  then  interact  with 
the  particle,  resulting  in  an  outgoing-scattered  wave  as  shown  in  figure  1.1.  The 
rigorous  solution  of  the  problem  requires  the  matching  of  boundary  conditions 
at  the  surface  of  the  particle  between  the  outside-incident  wave  and  the  inside 
wave,  and  the  inside  wave  and  the  outside-scattered  wave.  Therefore,  we  must 
investigate  the  nature  of  the  waves  in  the  medium  as  well  as  the  nature  of  the 
boundary  conditions.  We  do  this  in  the  next  three  subsections. 

I. 2.1  Maxwell’s  Equations  in  the  Scattering  Medium 

To  find  the  nature  of  the  electromagnetic  waves  in  the  medium,  we 
must  solve  Maxwell’s  equations  as  they  are  expressed  in  the  medium  itself. 
Consider,  then,  a  plane  wave  of  random  polarization,  and  a  spherical  particle 
with  which  the  wave  interacts,  as  in  figure  1.1.  In  the  scatterer  itself,  or  put 
differently,  in  a  macroscopic  medium  with  surface  charge  density  p  and  current 

J,  Maxwell’s  equations  are 


V  •  D  =  4  np 


(1.7) 
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1  BB 

v  X  E  =  — — 
c  at 

(1.8) 

V  •  B  =  0 

(1.9) 

1  dD  4x 

V  X  H  =  ~—  +  —  J 

c  at  c 

(1.10) 

where  H  is  the  magnetic  field,  B  is  the  magnetic  induction,  E  is  the  electric 
field,  and  D  is  the  electric  displacement.  The  electric  displacement  and  electric 
field  are  related  by 

D  =  eE,  (1.11) 

while  the  magnetic  field  and  magnetic  induction  are  given  by 

H  =  -B.  (1.12) 

Finally,  the  current  J  and  the  electric  field  E  are  related  by  Ohm’s  Law: 

J  =  aE,  (1.13) 

where  a  is  the  conductivity  of  the  medium.  In  section  1.2  we  assumed  that,  be¬ 
fore  interaction  with  the  scatterer,  the  magnetic  and  electric  fields  form  a  plane 
wave.  We  also  found  that  they  are  periodic  in  nature,  given  by  equations  1.5 
and  1.6.  When  they  interact  with  the  scatterer,  the  expression  for  the  E  and 
H  fields  changes  (as  we  shall  see).  Nonetheless,  they  maintain  their  periodic 
nature.  Thus,  in  the  scatterer,  the  magnetic  and  electric  fields  are  of  the  form 

E  =  E0(x)ex  p,u,t  (1.14) 

H  =  H0{x)expiut,  (1.15) 

where  E0{x)  and  H0(x)  are  the  complex  amplitudes.  By  taking  the  time 
derivatives  of  equation  1.14  and  equation  1.15  above,  plugging  into  the  two 
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curl  equations  (equations  1.8  and  1.10)  and  simplifying,  we  obtain 


V  X  E  =  —ikH 


(1.16) 


and 


V  X  H  =  ikm2E, 


(1.17) 


where  J  =  a E  and  D  =  eE  have  been  used.  The  new  term  m  in  equation  1.17 
is  called  the  complex  index  of  refraction ,  given  by2 


2  47 Tier 

m  =  e - , 

u 


(1.18) 


while  k  —  From  equation  1.18  it  is  apparent  that  the  complex  index  of 
refraction  is  a  function  of  frequency.  However,  in  practical  applications,  “m 
cannot  generally  be  determined  from  the  static  values  of  c  and  a  but  should 
be  determined  by  measurements  at  the  circular  frecuency  [van  de  Hulst  81, 
page  116]. 

The  curl  equations  remain  coupled.  However,  by  applying  the  di¬ 
vergence  equations,  we  can  decouple  them  as  follows.  From  equation  1.17  we 
get 

Taking  the  curl, 

VXE  =  ^VX(V><H)’ 


2In  so  defining  m,  we  have  separated  the  complex  part  of  the  problem  from  the  wave 
number  k.  Other  authors,  Jackson  in  particular,  prefer  to  include  the  complex  factor  in  the 
k  itself  [Jackson  75,  page  286].  Thus,  what  we  denote  as  mk  is  equivalent  to  Jackson’s  k. 
Furthermore,  our  expression  for  m  will  not  involve  the  parameter  /i,  since  we  shall  restrict 
our  analysis  to  substances  for  which  n  =  1. 
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or,  using  V  •  H  =  •  B  =  0, 

Vx£  =  -^-VJJf,  (1.19) 

Arm 2 

Equating  equation  1.16  to  equation  1.19  and  rearranging  terms,  we  get  an 
expression  involving  only  the  magnetic  field  H : 

V3H  -1-  lc2m2H  =  0.  (1.20) 

A  similar  (decoupled)  equation  can  be  obtained  for  the  electric  field,  E.  Thus, 
each  of  the  cartesian  coordinates  of  E  and  H  satisfy  the  scalar  wave  equation 

V2u  +  Ar2m2u  =  0.  (1.21) 

The  simplest  solution  to  the  scalar  wave  equation  is  very  enlightening.  This 

simple  solution  is 

u{z,t)  =  e-,fcmz+*wt. 

Notice  that,  if  m  has  an  imaginary  component,  then  the  wave  in  the  medium 
of  propagation  is  damped.  On  the  other  hand,  if  m  is  strictly  real,  the  wave 
suffers  no  attenuation,  although  it  does  undergo  a  phase  shift.  As  enlightening 
as  this  form  of  the  solution  might  be,  however,  it  does  not  represent  the  solution 
of  the  wave  in  a  spherical  particle,  nor  the  solution  of  the  wave  scattered  by  a 
small  sphere.  The  actual  solution  involves  spherical  waves,  which  is  the  topic 
of  the  next  section. 

1.2.2  Solution  of  the  Wave  Equation  in  Spherical  Coordinates 

Consider  the  spherical  coordinates  shown  in  figure  1.2.  Assuming  that 
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Figure  1.2:  A  spherical  coordinate  system.  The  origin  of  the  system  is  located 
at  the  center  of  the  scattering  sphere  [Arfken  85,  page  104]. 
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the  scalar  wave  equation  is  separable  in  r,  9  and  <j>,  its  solution  may  be  written 
as 

u(r,M)  =  «(r)0(0W). 

Using  the  laplacian  for  spherical  coordinates,  separating  the  radial  equation 
from  the  angular  dependence,  and  solving  the  resulting  three  equations,  we 
obtain3 

u(r,  6,  <{> )  =  ( A  cos (l<j>)  +  B  sin(/<^))P£(cos  9)zn(mkr ), 
where  A  and  B  are  constants,  /  and  n  are  integers  such  that 

n  >  /  >  0, 

P^(c os0)  is  the  associated  Legendre  polynomial,  and  zn(mkr)  is  the  spherical 
Bessel  function. 

Since  the  scalar  vector  function  represents  any  one  of  the  components 
of  the  vector  wave  function,  it  follows  that  the  solution  of  the  scalar  function 
must  be  related  to  the  solution  of  the  vector  wave  function.  Mie  shows  that,  if 
it  is  a  solution  of  the  scalar  wave  function,  then  the  vectors  Mu  and  Nu,  defined 

by 

Mu  =  Vx(ru)  (1.22) 

and 

mkNu  =  V  x  Afu,  (1.23) 

satisfy  the  vector  wave  equation.  Furthermore,  Mu  and  Nu  are  related  by 

mkMu  =  VxNu  (1.24) 


3The  separation  of  variables  in  spherical  coordinates  is  treated  in  many  textbooks.  See, 
for  instance,  [Slater  47,  pages  148  to  151]  and  [Arfken  85,  page  105  and  pages  448  to  450]. 


11 


Finally,  if  u  and  v  are  two  independent  solutions  to  the  scalar  wave  function, 
then  the  electric  and  magnetic  fields  are  given  by 

E  =  Mv  +  iN,  u  (1.25) 

and 

H  =  m(-Mu  +  iNv),  (1.26) 

respectively  [Mie  08,  page  120].  Thus,  given  two  independent  solutions  to  the 
scalar  wave  equation  of  a  wave  in  an  absorbing  medium,  we  can  find  the  electric 
and  magnetic  fields  which  satisfy  the  vector  wave  equation.  However,  the  wave 
equation  need  not  be  limited  to  a  dissipative  medium.  Equations  1.22  through 
1.26  apply  equally  well  to  an  incoming  ( incoming  from  the  point  of  view  of  a 
scatterer)  plane  wave  in  vacuum,  as  well  as  to  an  outgoing,  scattered,  spherical 
wave,  provided  m  is  appropriately  expressed.  The  exact  relationship  between 
these  waves  is  given  by  matching  the  appropriate  boundary  conditions  at  the 
skin  of  the  scatterer.  The  nature  of  these  boundary  conditions  is  tu"'  subject 
of  the  next  section. 

1.2.3  Boundary  Conditions  Between  Different  Media 

Equations  1.7  through  1.10  are  the  differential  form  of  Maxwell’s  equa¬ 
tions  for  waves  in  a  macroscopic  medium.  By  applying  the  divergence  and 
Stoke’s  theorem,  Maxwell’s  equations  can  be  written  in  their  integral  form: 

j>  D'nda  =  4tt  J  pcPx  (1.27) 

j>  B-nda  =  0  ( 1 .28) 
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D, 


Figure  1.3:  Infinitesimal  Gaussian  pillbox  used  to  determine  boundary  condi¬ 
tions  [Griffiths  81,  page  280]. 


1  dD 
-J+—7- 

c  c  at 


nda 


l  E-d\  =  --  /  ~nda , 
Jc  c  Js  at 


(1.29) 

(1.30) 


where  equations  1.27  and  1.28  are  defined  over  any  closed  surface  S  with  vol¬ 
ume  V,  and  equations  1.29  and  1.30  are  over  any  surface  S  bounded  by  the 
closed  loop  C.  In  this  form,  Maxwell’s  equations  are  particularly  useful  in  de¬ 
termining  the  boundary  conditions.4  Consider  a  very  thin  Gaussian  pillbox  as 
in  figure  1.3.  As  the  thickness  of  the  pillbox  goes  to  zero,  equation  1.27  leads 
to 

(£>2  -  DO-n  =  4?r£,  (1.31) 


4Griffiths  gives  a  particularly  good  presentation  of  the  limiting  process  used  to  arrive  at 
the  boundary  conditions  [Griffiths  81,  pages  280  and  281,  and,  on  a  different  context,  pages 
78  and  79] . 
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Figure  1.4:  Infinitesimal  Stokesian  loop  used  to  determine  boundary  conditions 
[Griffiths  81,  page  281]. 

where  E  is  the  surface  charge  density.  Meanwhile,  equation  1.28  leads  to 

(B2  -  Bx)  n  =  0.  (1.32) 

Similarly,  consider  the  Stokesian  loop  in  figure  1.4.  As  the  height  goes  to  zero, 
equation  1.30  leads  to 

nX(E2  —Ei)  =  0,  (1.33) 

while  equation  1 .29  leads  to 

nx(H2-H1)  =  —Kf ,  (1.34) 

c 

where  Kj  is  the  idealized  surface  current.  If  medium  1  has  an  index  of  re¬ 
fraction  m j  and  medium  2  has  an  index  of  refraction  m2,  both  finite,  then  the 
surface  current  Kj  vanishes  [van  de  Hulst  81,  page  117].  Then,  equation  1.34 
may  be  written  as 


nx{H2-H  0  =  0 


(1.35) 
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Equations  1.31,  1.32,  1.33  and  1.35  are  the  boundary  equations  to  which  the 
incoming,  inside  and  outgoing  wave  are  subject.  To  find  the  scattered  electric 
and  magnetic  fields  at  any  point  in  space,  all  that  remains  to  be  done  is  to  find 
the  waves  in  terms  of  known  quantities  and  apply  the  boundary  conditions. 
We  do  this  in  the  next  section. 

1.2.4  The  Mie  Coefficients:  Parameters  in  terms  of 
Known  Quantities 

A  review  of  the  progress  achieved  seems  appropriate  at  this  point5.  So 
far,  we  have  found  an  expression  for  the  plane  waves  in  vacuum  (equations  1.5 
and  1.6)  which  then  interact  with  the  spherical  scatterer.  We  have  found 
(equations  1.25  and  1.26)  the  expression  for  the  electric  and  magnetic  fields 
inside  and  outside  the  scatterer  in  terms  of  two  linearly  independent  scalar 
functions  u  and  v  (these  scalar  functions  are  to  be  determined).  We  have  also 
discussed  the  boundary  conditions  which  the  incoming  wave,  the  inside  and 
the  scattered  wave  must  satisfy  (equations  1.31,  1.32,  1.33  and  1.35).  In  this 
last  section  of  this  introductory  chapter,  we  apply  the  boundary  conditions 
to  obtain  the  solution  to  the  scattering  problem  in  terms  of  known  quantities. 
First,  however,  we  must  find  the  expression  for  the  scalar  functions  u  and  v.  As 
explained  in  [van  de  Hulst  81,  pages  121  through  123]  the  electric  and  magnetic 
fields  may  be  expressed  in  terms  of  the  scalar  functions  u  and  v  as  follows. 

5The  credit  goes  to  Mie,  of  course.  In  this  section  we  continue  to  summarize  his  theory, 
as  presented  by  van  de  Hulst. 
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The  outside,  incident  wave: 

«  =  «'"‘cos((*)f;(-ir-rLiTTf,'(cos«)j„(tr),  (1.36) 

n=x  n(n  +  I) 

OO  <y  I  1 

v  =  e"‘jin(*)  £(-i)*— -1—  P>(cos«)j,(ir).  (1.37) 

n=i  n(n  +  l) 

The  angles  <£  and  #  are  as  defined  in  figure  1.2,  n  is  an  integer  from  0  to  oo,  jn 
is  the  spherical  Bessel  function  of  the  first  kind  and  P,j(cos0)  is  the  associated 
Legendre  polynomial. 

The  inside  wave: 

u  =  e,wtcos (<{>)  rncn(-i)n  ^  *  *  Pjfcosfl^mfcr),  (1.38) 

n=i  n(n  +  1) 

V  =  e'^sin (<t>)^2mdn(-i)n  ^  *  j  P,,1  ( cos0)jn (mkr),  (1.39) 

n=i  n(n  +  l) 

and  the  outside,  scattered  wave 

OO  ey  i  i 

u  =  eiu,t cos(<j>)  ^2  —an(—i)n—  ^——P*(cos9)h<£\kr),  (1.40) 

n=i  n(n  +  l) 

v  =  e,w*sin(<£)  ^  — 6n(— i)n-^-i~P^(cos^)/i^(^r),  (1.41) 

„=i  n(n  +  l) 

where  h^(kr)  is  the  spherical  Hankel  function  of  the  second  kind.  In  the 
equations  for  the  inside  and  the  scattered  waves,  (equations  1.38  through  1.41) 
Am  &n>  Cn  and  dn  are  called  the  Mie  coefficients.  These  coefficients  are  the 
sought-after  parameters  which  are  expressed  in  terms  of  known  (or  functions 
of  known)  quantities.  By  matching  the  appropriate  boundary  conditions  at  the 
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skin  of  the  particle  (see  section  1.2.3),  we  obtain  for  the  Mie  coefficients 

Sn(y)Sn(*)  ~  mSn(y)S'n(x) 


an  = 


bn  = 


S'n(y)(n(x)  ~  mSn{y)Cnix)  ’ 

™S'n{y)Sn(x)  -  Sn(y)S'n(x) 
mS'n(y) C„(x)  -  Sn(y)Cn(x)  ’ 


(1.42) 


(1.43) 


where  x  —  ka  —  y  =  mka,  and  m  is  the  complex  index  of  refraction  of  the 
medium.  5„  and  Cn  are  the  Riccati- Bessel  functions  defined  by 


Sn(z)  =  zjn(z),  (1-44) 

and 

Cn(z)  = -znn(z).  (1.45) 

where  jn(z)  and  nn(z)  are  the  spherical  Bessel  functions  of  the  first  and  second 
kind,  respectively6.  The  function  ^(n)  is  defined  in  terms  of  S„  and  Cn: 

C(z)  =  Sn(z)  +  iCn(z). 

Thus,  for  a  given  radius  r,  a  wavelength  A,  and  a  complex  index  of  refraction 
m,  we  can  find  the  corresponding  Mie  coefficients.  Having  found  the  Mie 
coefficients,  we  can  then  find  the  fields  at  any  point.  However,  the  distribution 
of  the  size  of  the  particles  in  the  cloud,  the  shape  of  the  cloud  and,  most 
importantly,  the  direction  in  which  the  detector  is  “looking,”  make  the  problem 
nontrivial.  In  the  next  chapter  we  investigate  these  factors. 


6See  section  3.1  for  the  generating  function  used  for  the  Riccati- Bessel  functions. 


Chapter  2 


Particle  Distributions  and  Cloud  Shapes 

2.1  Scattered  and  Absorbed  Signals  in  Clouds  with 
Nonuniform  Particle  Distributions 

In  problems  of  practical  importance,  the  signal  is  affected  by  a  cloud 
of  particles,  as  opposed  to  just  one  particle.  Recall  the  problem  at  hand:  a 
signal  may  be  masked  by  the  light  from  a  strong  source  scattered  into  the 
detector,  or  the  signal  itself  may  be  attenuated  by  the  cloud.  Although  the 
cloud  may  have  any  geometrical  shape,  the  particles  tend  to  be  spherical.1  If 
we  further  assume,  for  now,  that  all  particles  in  the  cloud  have  radius  a  and 
that  the  number  of  particles  per  unit  volume  is  N,  then  the  light  intensity 
scattered  per  unit  volume  in  the  direction  9,  is 


Itotai(0,a)  =  NI(9,a) 


(2.1) 


where  I{9,a)  is  the  intensity  of  light  scattered  by  one  particle.  In  terms  of  the 
Mie  coefficients,  the  intensity  is  given  by2 


/(M)  =  /0 


+  *2(0) 

2  fc2r2 


(2.2) 


'This  certainly  seems  true  of  particles  produced  by  rocket  motors.  See,  for  instance, 
[Dawbarn  80] . 

2On  first  inspection,  it  seems  as  though  the  intensity  1(9,  a)  is  not  a  function  of  the  radius 
a  of  the  particle.  However,  we  must  keep  in  mind  that  the  Mie  coefficients,  which  are  a  part 
of  »i  and  »2  by  equations  2.2  through  2.6,  are  radius-dependent. 
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where  r 

is  the  distance  from  the  detector  to  the  scatterer,  /„  is 

the  intensity 

of  the  incident  wave,  and  tj(0)  and  12(0)  are  the  scattering  functions  for  the 
perpendicular  and  plane  polarized  scattered  waves,  respectively.  They  are  given 

by 

*‘i(«)  =  |Si(«)|2, 

(2.3) 

and 

;2(*)  =  |s2(0)|2. 

(2.4) 

In  turn,  S\(0)  and  «S2($)  are  given  by3 

Sl(0)  -  X)  /  ,  ,x{anMcos0)  +  6nrn(cos0)}, 

n=l  n(n  +  !) 

(2.5) 

and 

OO  rt  |  I 

Si(e)  =  J2  n(n  +  ^{^(costf)  +  anr„(cos0)}. 

(2.6) 

Finally, 

7rn(cos0)  and  rn(cos  9)  are  given  by 

/  /.x  dPn(cos9) 

x„cos0)-  , 

a  cost) 

(2.T) 

and 

rn(cos0)  =  cos(0)7rn(cos  9)  —  sin2  , 

a  cos  9 

(2.8) 

where  Pn(cos9)  is  the  Legendre  polynomial  of  order  n.  The  second  part  of  the 
problem  involves  the  attenuation  of  the  signal  itself.  This  attenuation  is  given 

by 

I  =  I0e~'R,  (2.9) 

3The  functions  S\  and  S 2  are  not  to  be  confused  with  the  Riccati- Bessel  functions  S„ 
defined  by  equation  1.44. 
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where  R  is  the  distance  the  signal  has  traveled  in  the  absorbing  cloud,  and  7  is 
the  extinction  coefficient ,  which  is  a  function  of  the  size  of  the  particles,  a,  the 
number  of  particles  per  unit  volume,  N,  and  the  extinction  efficiency  factor , 

Q ext- 

7  =  Nna2Qext.  (2.10) 


In  terms  of  the  Mie  coefficients,  the  extinction  efficiency  factor  is  given  by 


Q'xt  =  £(2n  +  l)Re(a„  +  6„) 


(2.11) 


Finally,  in  the  above  equation,  x  is  the  ratio  of  the  particle’s  circumference  to 
the  wavelength  of  the  light: 

*  =  ^  (2.12) 

Realistically,  however,  the  particles  in  the  cloud  will  be  of  different 
sizes.  The  number  of  particles  per  unit  volume  of  a  given  radius  is  determined 
by  a  distribution  function.  The  number  density  of  particles  with  radii  between 
a  and  a  +  da  (as  given  by  the  distribution)  replaces  N  in  formulas  2.1  and  2.10. 
The  effect  of  the  whole  cloud  is  obtained  by  adding  the  contributions  of  all  the 
individual  parts.  For  instance,  for  the  total  intensity  of  the  light  scattered  in 
the  direction  0,  the  expression  is 


W0)  =  E'.(M,) 


(2.13) 


where  7,(0,  a,)  is  the  intensity  of  light  scattered  at  an  angle  9  by  particles  of 
radius  a j.  Similarly,  the  total  attenuation  coefficient  7  becomes 


7  =  /  Q'M^(a)*a2da, 

JO 


(2.14) 
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where  ^ (a)  represents  the  normalized  particle  size  distribution: 


a)da ,  then,  represents  the  number  of  particles  with  radii  between  a  and 
a  +  da  per  unit  volume.  In  the  next  section  we  will  find  the  actual  particle 
population  for  a  given  distribution. 

2.2  Determination  of  Particle  Population 

Throughout  this  paper,  we  have  summarized  the  results  of  the  Mie 
theory  and  treated  the  scattering  of  an  electromagnetic  signal  in  very  general 
terms.  We  will  now  specialize  the  discussion  to  the  cloud  of  aluminum  oxide 
particles  found  in  the  exhaust  of  rockets.  The  aluminum  oxide  is  produced  as 
a  byproduct  of  aluminized  solid  fuel  used  to  propel  the  rockets.  As  the  fuel 
burns,  the  aluminum  particles  in  the  propellant  vaporize  and  react  with  oxygen 
to  form  molten  aluminum  oxide  [Mularz  72].  This  example  is  of  current  interest, 
for  it  applies,  among  other  problems,  to  the  recovery  of  satellites  by  the  Space 
Shuttle.  One  of  the  missions  of  the  Shuttle  is  to  recover  satellites  before  they  fall 
back  to  the  Earth.  If  the  satellite  has  not  been  visually  detected  (or  cannot  be 
visually  detected  due  to  distance,  operational  constraints,  etc.),  then  it  must  be 
detected  by  the  infrared  signature  the  satellite  possesses.  If  the  Shuttle  is  firing 
its  rockets  to  maneuver,  then  the  signal  may  have  to  pass  through  the  cone- 
shaped  cloud  of  the  Shuttle’s  exhaust.  The  size  of  the  particles  in  the  cloud 
will  not  be  of  one,  and  only  one,  radius.  They  will,  instead,  be  distributed 
throughout  a  wide  range  of  values.  However,  to  illustrate  the  method  used  to 
obtain  the  solution  for  a  non-constant  distribution,  we  first  look  at  the  simplest 
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of  distributions:  a  constant  distribution.  Assuming  there  are  M  particles  of 
density  p  and  constant  radius  a ,  the  total  mass  M  of  the  particles  is4 


M  =  Af-ira^p. 


Solving  for  N,  the  total  number  of  particles, 


M 


§7ra3p 


This  result  is  most  useful  if  expressed  as  a  particle  density: 


A r  =  _M _ 1 _ , 

*7ra3/)  V olumedoud' 

where  N  represents  the  number  of  particles  per  unit  volume  and  V olumeciou<i 
depends  on  the  shape  of  the  cloud. 


For  a  non-constant  distribution,  we  must  take  a  weighted  average. 
Consider,  for  example,  the  following  distribution,  which  is  the  distribution  of 
particles  produced  by  the  Titan-IIIC  rocket  exhaust5: 


^Jf(a)  =  0.012(2a)2  exp[-1.89  x  10"2(2a)2],  (2.15) 


which  is  shown  in  figure  2.1.  By  definition  of  a  distribution,  AA/*(a,)  is  the 
number  of  particles  with  radii  between  a,  and  a*  +  Aa.  Let  uo/(a,)  be  the 


^Notice  that  the  actual  numierof  particles  is  denoted  by  M,  while  the  number  density  is 
denoted  by  N. 

5It  may  seem  as  though  we  Me  using  one  system  (the  Titan-IIIC  rocket)  to  talk  about 
another  system  (the  Space  Shuttle).  However,  the  reader  must  bear  in  mind  that  the  Shuttle 
was  introduced  only  as  an  example.  Certainly,  there  must  be  other  examples  to  which  this 
analysis  applies.  Furthermore,  the  distribution  of  the  Shuttle’s  exhaust  must  resemble,  at 
least  in  shape,  the  distribution  presented  in  figure  2.1 
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Figure  2.1:  Particle  distribution  for  a  Titan-IIlC  rocket  motor  [Dawbarn  80, 
page  126]. 


volume  of  the  particles  with  radii  between  a,  and  a,  +  A  a.  Then  the  mass  of 
all  the  particles  present  is 


M  =  Y'  vol(a,i  )pAjV  (a, ) 

i 

M  =  pT>o/(q,-)A^(at). 

t 

Solving  for  AA/’(aI),  we  get 


A^(a<)  = 


M 


(2.16) 


p52iVol(ai)' 

We  have  assumed,  for  simplicity,  that  all  the  burned  mass  turns  to  smoke 
particles  or  other  absorbing  particles.  This  is  usually  not  the  case.  Instead, 
only  a  fraction  /  of  the  amount  burned  turns  into  particles.  Thus  equation  2.16, 
expressed  as  a  number  density,  becomes 

/M  1 


A  N(ai)  = 


(2.17) 


P  Ei  vol(ai)  Volumeciou<i  ’ 
where  A/V(a,)  is  the  number  of  particles  per  unit  volume  with  radii  between 
a,  and  a;  +  Aa.  To  encode  the  above  ideas  into  a  computer  program,  we  need 
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to  write  a  program  to  solve  the  scattering/ absorption  problem  for  a  “constant” 
distribution.  We  then  take  care  of  the  changing  nature  of  the  distribution  by 
incrementing  the  radius  a  small  amount,  and  solving  that  “constant”  distribu¬ 
tion  problem.  However,  the  amount  of  light  (scattered  or  absorbed)  reaching 
the  detector  depends  not  only  on  the  distribution  of  particles  in  the  cloud,  but 
also  on  the  shape  of  the  cloud  itself.  We  investigate  the  effect  of  the  cloud  on 
the  signal  in  the  next  section. 

2.3  Effect  of  the  Cloud  Shape 

We  now  turn  to  the  question  of  the  effect  the  shape  of  the  cloud  has 
on  the  scattering  and  absorption  of  light.  Again,  we  start  with  a  very  simple 
model:  a  detector  at  the  center  of  a  spherical  cloud  composed  of  constant 
size  particles.  Figure  2.2  shows  the  arrangement.  In  this  case  there  is  perfect 
symmetry.  The  intensity  of  the  scattered  light  reaching  the  detector  is  the 
intensity  of  the  light  in  the  solid  angle  subtented  by  the  cone: 

/flc  _  rh  i(0) 

I  ~  2tt  /  r2dr  I  sin  Od0-±f-.  (2.18) 

Jo  JQ\  r 

Substituting  equation  2.2  for  the  intensity  i(0)  in  equation  2.18,  we  can  express 
the  ratio  of  the  scattered  intensity  to  the  original  intensity  as 

{  =  HF5  C 8in(#)(il  m  +  (2' 19) 

where  0  is  the  angle  the  detector  makes  with  the  z  axis,  0'  is  the  scattering 
angle  ,  0'  =  n  —  0,  and  A0  =  02  —  0\  is  the  resolution  angle  of  the  detector 
(see  figure  2.2).  Recall  that  i’i(0)  and  i2(0)  are  functions  of  the  Mie  coeficients 
(see  equations  2.3  through  2.6).  Thus,  all  the  parameters  in  equation  2.19  are 
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t  tt  t  t 

Figure  2.2:  A  detector  in  a  spherical  cloud.  The  intensity  reaching  the  detector 
is  the  intensity  of  light  in  the  solid  angle  subtended  by  the  detector’s  field  of 
view. 
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known.  Consequently,  for  a  given  size  particle,  we  now  have  an  expression 
for  the  scattered  intensity  (the  absorption  of  the  signal  remains  as  given  by 
equations  2.9  through  2.12).  For  a  distribution  of  particles  with  n  different 
sizes,  we  simply  repeat  the  process  n  times,  each  time  with  the  appropriate 
particle  radius. 

For  a  cone-shaped  cloud,  the  problem  is,  in  many  respects,  similar 
to  the  problem  of  the  spherical  cloud — and  uniquely  different  in  others.  For 
instance,  if  the  solid  angle  subtented  by  the  detector’s  view  angle  is  com¬ 
pletely  contained  within  the  conical  cloud,  then  the  problem  is  identical  to 
the  spherical-cloud  problem-with  the  provision  that  the  particle  density  will 
be  higher  for  the  cone  cloud  because  of  its  smaller  volume  for  compatible  di¬ 
mensions  (see  figure  2.3).  However,  if  the  detector  is  “looking”  off  the  axis 
of  the  cloud,  as  in  figure  2.4,  then  the  problem  looses  all  resemblance  to  the 
spherical  cloud  problem  and  becomes  much  more  complicated.  Specifically,  the 
intensity  of  the  scattered  light  will  depend  on  the  following  factors: 

1.  The  relative  sizes  of  the  solid  angle  of  the  detector’s  field  of  view  and  the 
conical  cloud. 

2.  The  angle  the  axis  of  the  cloud  makes  with  the  sun’s  rays  (the  z  axis). 

3.  The  viewing  angle  of  the  detector  with  respect  to  the  axis  of  the  cloud. 
This  angle  does  not  have  to  be  in  the  same  plane  as  the  sun’s  rays  and 
the  axis  of  the  cloud. 

We  treat  the  case  in  which  conditions  1  and  2  above  apply,  with  the  detector’s 
view  axis  on  the  z-y  plane,  as  in  figure  2.3.  A  detector  view  width  of  two 
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Figure  2.3:  The  cloud  as  a  cone.  The  angle  subtended  by  the  detector  is  smaller 
than  the  angle  subtended  by  the  cloud 
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Figure  2.4:  The  detector’s  view  axis  off  the  y-z  plane. 

degrees  is  used,  while  the  ratio  of  the  length  of  the  cone  cloud  to  the  radius 
of  the  base  is  set  as  10  to  1,  giving  a  cloud  angle  of  11  degrees.  Thus,  the 
detector’s  view  angle  is  contained  within  the  cloud. 


Chapter  3 


Data  Used,  Results  and  Calculations 

3.1  Data  and  Method  of  Calculation 

The  computer  program  can  be  described  as  a  series  of  iterations  of  the 
equations  in  chapters  1  and  2.  Specifically,  it  divides  the  distribution  function 
of  the  particles  (smoke  or  aluminum  oxide)  into  50  equally  spaced  regions.1 
Thus,  since  it  is  assumed  that  the  smallest  radius  in  the  distribution  is  zero 
and  that  the  largest  radius  is  12  /im,  the  center  of  the  zth  division  corresponds  to 
a  radius  of  a,  =  (i  +  -5)(0.24)  /zm,  or  to  a  diameter  of  D,  =  (z  +  ^)(0.48)  pm.  For 
each  radius  thus  calculated,  it  finds  the  Mie  coefficients  using  equations  1.42 
and  1.43.  The  Riccati- Bessel  functions  used  in  equations  1.42  and  1.43  are 
generated  using  the  following  expressions,  taken  from  [Gumprecht  51,  pages 
xii  and  xiii] : 


SW(X)  = 

— — i-S.-i(i)  -  s„-2(i), 

X 

(3.1) 

2n  —  1  .  .  .  . 

Cn(x)  = 

- -Cn.tix)  -  Cn-2(X), 

X 

(3.2) 

given  that  the  two  initial  functions  are  given  by 

S0(x)  =  sin(x),  Si(x)  =  -  cos(x), 


1  We  found  that  50  iterations  give  results  accurate  to  one,  and  for  some  wavelengths,  two 
significant  figures.  Since  the  results  are  reported  as  graphs,  we  felt  that  any  higher  accuracy 
would  be  superfluous. 
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and 

Co(x)  =  cos(x),  C\(x)  =  co8jJ^  +  sin(x). 

The  derivatives  of  the  Riccati- Bessel  functions  follow  from  the  above  definitions. 
Having  found  the  Mie  coefficients,  the  program  then  calculates  the  scattering 
functions  for  the  perpendicular  and  plane  polarized  scattered  waves,  i\(9')  and 
i2(0')»  as  defined  by  equations  2.3  through  2.8.  Recall  (section  2.3)  that  the 
scattering  angle  9'  and  the  angle  the  detector  makes  with  the  z  axis  (the  sun 
rays)  9  are  related  by 

9'  =  x  -  9. 


The  program  then  calculates,  by  use  of  equation  2.19,  the  ratio  of  the  intensity 
of  the  scattered  signal  to  the  intensity  of  the  incoming  signal,  -f-,  for  the  view 
angle  9  from  0  to  180  degrees.  This  corresponds  to  a  scattering  angle  of  180  to 
0  degrees.  The  program  repeats  the  above  process  for  all  the  divisions  in  the 
distribution. 


In  calculating  the  attenuation,  we  assumed  that  the  signal  travels  the 
extent  of  the  cloud.  Thus,  the  attenuation  has  no  angular  dependence.  How¬ 
ever,  the  contributions  of  all  the  partitions  must  still  be  added,  as  indicated  by 
equation  2.14.  In  addition,  we  found  that  the  range  of  values  in  the  attenua¬ 
tion,  as  given  by  equation  2.14,  is  very  wide.  Thus,  we  express  it  in  terms  of 
decibels: 


attenuation  =  10 


log 


h 

r 


or 


attenuation  =  4.3437  Rc. 
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Figure  3.1:  The  imaginary  part  of  the  index  of  refraction  for  aluminum  oxide 
at  2950  K  [Bakhir  77]. 

To  calculate  the  Mie  coefficients,  we  need  the  value  of  the  index  of 
refraction  of  the  substance  at  the  frecuency  (and,  thus,  the  wavelength)  in 
question.  We  used  the  values  calculated  by  Bakhir  et  al.  [Bakhir  77],  who 
published  a  graph  of  the  dependence  of  the  imaginary  part  of  the  index  of 
refraction  for  aluminum  oxide  at  2950  K  as  a  function  of  wavelength.  We  have 
reproduced  their  results  in  figure  3.1.  They  also  published  the  dependence  of 
the  real  part  of  the  index  of  refraction  on  the  wavelength  in  the  form  of  a  table. 
We  reproduce  this  result  in  figure  3.2.  For  the  index  of  refraction  of  smoke,  we 
used  m  =  1.75  —  i0.45,  which  was  taken  from  Ramaswamy  [Ramaswamy  85]. 

It  was  assumed  that  1  kg  of  fuel  was  burned,  and  that  the  fraction  of 
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mass  burned  that  turned  into  aluminum  oxide  or  smoke  particles  is  0.05.  The 
densities  were  taken  as  2.0  x  103  kg/m3  for  smoke  and  as  2.7  x  103  kg/m3  for 
aluminum  oxide.  The  results  are  shown  in  the  next  section. 

3.2  Results  of  the  Models 

The  results  are  divided  into  two  groups  of  graphs.  The  first  group, 
figures  3.3  through  3.6,  show  the  intensity  (expressed  as  the  ratio  of  the  in¬ 
tensity  reaching  the  detector  and  the  initial  intensity,  j-)  of  scattered  light  the 
detector  receives.  Each  one  of  the  graphs  show  the  intensity  as  a  function  of 
the  cloud  extention  ( Rc  in  figure  2.3)  and  the  detector’s  view  angle  for  a  given 
wavelength  and  index  of  refraction.  For  instance,  figure  3.3  represents  the  ratio 
for  the  wavelength  A  =  4.04  /«m  and  index  of  refraction  m  =  1.76  —  z'0.0102 
as  a  function  of  viewing  angle  for  cloud  extentions  Rc  of  1,  10,  100  and  1000  me¬ 
ters.  The  scattering  particles  are  aluminum  oxide  in  figures  3.3  through  3.5  and 
smoke  in  figure  3.6.  It  was  assumed  that  the  size  distribution  for  the  smoke 
particles  was  identical  to  the  distribution  of  aluminum  oxide.  The  informa¬ 
tion  presented  in  these  graphs  may  be  used  to  determine  the  extent  to  which 
the  scattered  light  masks  the  signal.  Ultimately,  the  ability  of  the  detector  to 
distinguish  the  signal  from  the  “noise”  introduced  by  the  scattered  light  will 
depend  on  the  intensity  of  the  incident  light,  /„,  and  the  intensity  of  the  signal 
itself,  I  signal'  While  the  intensity  of  sunlight  is  well  known,  a  determination  of 
the  signal  intensity  and  the  detector’s  ability  to  distinguish  the  signal  from  the 
scattered  “noise”  is  beyond  the  scope  of  this  project. 

The  second  group  of  results  is  shown  in  figure  3.7.  This  graph  shows 
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the  attenuation  (as  defined  in  section  3.1)  of  the  signal  itself  as  a  function  of 
extention  Rc  of  the  cloud.  Since  the  attenuation  in  aluminum  oxide  for  signals  of 
wavelengths  4.04  /x m,  6.25  /xm  and  8.0  /xm  was  found  to  be  essentially  identical, 
only  the  graph  for  A  =  8.0  /xm  is  shown.  In  addition,  we  also  show,  in  the  same 
graph,  the  attenuation  of  light  at  A  =  10.0  /xm  for  smoke.  From  this  graph  we 
see  that  the  attenuation  is  in  the  order  of  ldb  for  a  cloud  extention  of  10  meters, 
falling  rapidly  for  greater  cloud  demensions.  Thus,  attenuation  should  not  be 
a  factor  in  the  detection  of  signals  for  cloud  extentions  greater  than  10  meters. 
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3.3  Conclusions 

The  results  of  the  model  for  the  scattered  intensity,  figures  3.3  through 
3.6,  show  the  ratio  of  scattered  intensity  to  incident  intensity  •£.  This  infor¬ 
mation  can  be  used,  in  conjunction  with  information  on  the  intensity  of  the 
sunlight  and  the  signal  itself,  to  determine  the  effect  of  the  scattered  light  on 
the  detector’s  ability  to  detect  the  signal.  To  calculate  the  graphs,  we  assumed 
that  the  amount  of  fuel  burned  is  1  kg,  that  the  fraction  of  particles  converted 
to  aluminum  oxide  or  smoke  is  0.05,  and  that  the  densities  are  2.7  x  103  kg/m3 
for  aluminum  oxide  and  2.0  x  103  kg/m3  for  smoke. 

The  results  for  the  attenuation  are  much  easier  to  interpret.  Figure  3.7 
shows  that  the  attenuation  of  a  signal  in  the  range  of  10  /im  is  of  the  order  of 
1  db  for  a  cloud  of  10  m.  The  attenuation  decreases  rapidly  with  increasing 
cloud  dimensions.  Thus,  the  attenuation  should  not  be  a  problem  in  signal 
detection  once  the  cloud  has  dissipated  to  10  meters  or  greater. 

3.4  Recommendations  for  Further  Work 

Future  investigations  of  this  problem  should  look  at  the  problem  of 
a  detector  searching  for  a  signal  in  a  direction  not  in  the  plane  defined  by  the 
cloud  axis  and  the  sun’s  rays.  Although  the  theory  of  the  problem  is  the  same 
as  in  this  presentation,  the  unusual  geometry  of  an  off-plane  detector  makes  it 
much  more  complicated. 


Figure  3.3:  j-  for  the  aluminum  oxide  particle  distribution  of  a  Titan-IIIC 
rocket  at  A  =  4.04  fim. 


detector’s  view  angle  9  (degrees) 
(scattering  angle  9‘  =  7r  —  9) 


Figure  3.4:  -j-  for  the  aluminum  oxide  particle  distribution  of  a  Titan-IIIC 
rocket  at  A  =  6.25  /jm. 
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Figure  3.5:  j-  for  the  aluminum  oxide  particle  distribution  of  a  Titan-IIIC 
rocket  at  A  =  8.0  um. 
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Figure  3.6:  y  for  the  smoke  particle  distribution  of  a  Titan-IIIC  rocket  at 


Figure  3.7:  The  attenuation  for  the  aluminum  oxide  and  smoke  particle  distri 
bution  of  a  Titan- IIIC  rocket  at  A  =  8.0  n m  and  A  =  10.0  /im,  respectively. 


Appendix  A 


The  Computer  Program 


We  have  tried  to  make  the  program  as  self-explanatory  and  user- 
friendly  as  possible.  Experience  shows,  however,  that  what  is  clear  to  the 
programmer  need  not  be  clear  to  the  user.  Thus,  we  include  the  following 
“pointers”  in  the  hopes  of  clearifying  the  workings  of  the  program. 

•  Whenever  possible,  the  equations  or  formulas  in  the  program  are  written 
in  the  form  they  are  written  in  the  text.  Unfortunatelly,  this  is  not  always 
possible  due  to  the  limited  selection  of  characters  available  to  the  pro¬ 
grammer.  For  instance,  while  the  text  can  easily  handle  characters  such 
as  A  and  S\,  the  programmer  does  not  have  access  to  those  characters. 
In  those  cases,  we  either  spell  the  name  of  the  variable  (“lambda”  for  A, 
for  instance)  or  we  try  to  associate  a  discriptive  name  with  the  variable 
(“numdivisions”  for  number  of  divisions.) 

•  The  program  asks  the  user  to  input  the  wavelenght  of  the  radiation,  the 
number  of  partitions  in  the  distribution,  the  extention  of  the  cloud  and  the 
index  of  refraction.  When  supplying  the  data,  it  is  extremely  important 
that  the  data  be  in  the  units  specified  in  the  program;  otherwise,  although 
the  program  may  run,  the  results  will  be  meaningless  (the  program,  as 
stands,  does  not  check  for  the  correct  units  nor  dimensions.) 
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•  In  subroutines,  the  input  and  output  arguments  are  separated  in  one 
of  two  fashions:  by  a  space  between  the  last  input  argument  and  the 
first  output  argument,  or  by  placing  the  output  arguments  in  a  new  line. 
Additionally,  at  the  beginning  of  every  subroutine,  we  say,  in  a  comment 
box,  which  arguments  are  inputs  and  which  are  outputs. 

•  In  addition  to  output  arguments,  subroutines  may  also  create  data  files. 
These  data  files  are  then  used  by  other  subroutines. 

•  The  results  for  the  scattering  problem  (j^)  are  stored  in  the  file  sums.dat, 
while  the  attenuation  is  stored  in  the  file  data  file  attenuation.dat. 

Turn  to  the  next  page  for  the  program. 
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program  scattermain 

****♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 
♦Comment:  Throughout  the  program,  comments  will  be  enclosed 
*in  boxes  (such  as  this) .  This  program  asks  the  user  to 
♦input  the  wavelength.  For  the  program  to  run  correctly, 

*it  should  be  expressed  in  micrometers  (IE-06) .  On  the  other 
♦hand,  the  cloud  radius  should  be  in  meters.  It's  very  important 
*to  input  the  complex  index  of  refraction  as  follows:  if  z  is 
♦a  complex  number,  z  =  a  +  ib,  where  a  is  the  real  part  of  z  and 
*b  is  the  imaginary,  then  it  should  be  typed  as  (a,b) .  Thus, 

*z  *  1.75  -  i0.45  would  be  typed  as  (1.75,-0.45). 

♦ 

♦NOTE:  if  data  for  a  density,  fraction  of  burned  mass,  or  any 
♦other  parameter,  different  from  the  ones  used  in  the  program, 
♦is  wanted,  then  they  must  be  changed  in  the  declarations  below 
**********♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦*♦♦♦♦♦♦♦♦♦♦♦♦*♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

double  precision  left , right .step, numberdensity(0 : 1000) 
double  precision  Q(0 : 1000) ,totalgamma,gamma(0 : 1000) 
double  precision  lambda, viewwidth.cloudradius, PI 
double  precision  volume, k, intensity, sum(20) , a, n 
double  precision  density .massburned .cloudbase 
complex  refraction 
integer  i,numdivisions 


external  qamp , intensitysub 
external  particlepopulation 
parameter  (PI  *  3.141592635) 
parameter  (viewwidth  ■  2.0) 
parameter  (massburned  =  0.05) 
parameter  (left  =  0.0) 
parameter  (right  *  12.0) 

open  (l,file*'viewwidth' ,status*’new' ) 
write  (1 ,*) viewwidth 
close  (1) 

print  *, 'input  density  in  kg/m(3)' 
read  *, density 

print  *, 'input  number  of  divisions  (even)' 
read  * .numdivisions 

print  input  wavelenght  in  micrometers  (10  (-6)m) ' 
read  *, lambda 

print  *, 'input  cloudradius  in  meters' 
read  * .cloudradius 

print  *, 'input  the  complex  index  of  refraction 
in  this  format:  m  *  a+ib  *  (a.b)' 
read  *, refract ion 


k  *  (2 .0*PI)/(lambda*lD-06) 
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************************************************************** 
♦the  volume  is  different  for  different  cloud  shapes!!  If  data 
♦for  a  cloud  not  in  a  cone  shape  is  wanted  (say,  a  sphere), 
♦then  the  following  formula  must  be  changed 
************************************************************** 

cloudbase  =  cloudradius/10 . 0 

volume  ■  (1 .0/3 . 0)*PI*(cloudradius**2)*cloudbase 
totalgamma  =  0.0 


step  *  (right  -  left)/numdivisions 

call  part iclepopulat ion (numdivisions , left , step .volume , 
*  massburned, density, PI) 

open  (1,  file  *  'population',  status  =  'old') 
do  10,  i  *  0,  (numdivisions  -1),1 
read(l,*)  numberdensity (i) 

10  continue 

close(l,  status  =  'delete') 

do  09,  i  =  0,  (numdivisions  -1),  1 
n  »  (2*i  +  l)/2.0 
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a  *  left  +  (n  *  step) 

cadi  qamp (a/lambda, refraction,  Q(i)) 

gamma(i)  *  numberdensity(i)*PI*((a*lD-06)**2)*Q(i) 

totalgamma  *  totalgamma  +  gamma(i) 

call  intensity sub (k .numberdensity , cloudradius) 


open  (1,  file  *  ’intensity' ,  status  *  ’old’) 
do  11,  j-1,  19 

read  (1 ,*) intensity 
sum(j)  *  sum(j)  +  intensity 
11  continue 

close  (1,  status  *  ’delete’) 

print*, 'this  is  division  #  /  out  of  numdivisions : ’ 
write(*,*)i,numdivisions 
09  continue 

open(l,  file  *  ’totalgamma’,  status  *  'new') 
vrite(l ,*)totalgamma 
close(l) 

* 

*  comment:  if  needed,  the  total  coefficient  of  extinction 

*  (totalgamma)  may  be  obtained  from  the  file  totalgamma.dat 

*  above. 

* 
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open  (2,  file  *  'sums',  status  *  'new') 
do  13,  i  *  1,  19 
write  (2,*)sum(i) 

13  continue 

close  (2) 

open  (1,  file  *  'attenuation' .status  *  'new') 
write  (1,*) (4.343*totalgamma*cloudradius) 
close(l) 
end 

************************************************************* 
♦Particlepopulation  finds  the  actual  number  of  particles 
♦as  dictated  by  the  distribution  and  the  amount  of  fuel 
♦burned.  It  calls  on  the  function  "distribution." 

♦All  the  arguments  are  inputs.  The  output  is  the 
♦data  file  population.dat,  which  contains  the  normalized 
♦distribution  of  particles. 

************************************************************* 

subroutine  part iclepopulat ion (numdivisions .left .step, 
♦  cloudvolume , massburned , density , PI ) 


integer  numdivisions , i 


47 


double  precision  n,a,left .step, numnormalized(0: 1000) 
double  precision  particlevolume, denominator ,numatrmax 
double  precision  distribution, massbumed.totalnum 
double  precision  cloudvolume,numberdensity(0 : 1000) ,PI 
double  precision  sumofparticles 
double  precision  number(0 : 1000) 
external  distribution 

sumofparticles  *  0.0 
do  9,  i=0,(numdivisions-l),l 
n  *  (2*i  +l)/2.0 
a  *  left  +  (n*step) 
number (i)  *  distribution (a, step) 
sumofparticles  *  sumofparticles  +  number(i) 
write(* ,  *)  i , sumofparticles 
9  continue 


denominator  *  0.0 
do  10,  i*0, (numdivisions-1) ,1 
n  *  (2*i  +l)/2.0 
a  *  left  +  (n*step) 

numnormalized(i)  *  distribution(a,step)/sumofparticles 
particlevolume  *  4. 0/3.0  *  PI  *  (a* lD-06)**3 
denominator  *  denominator  +  numnormalized(i) 
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*  *particlevolume 

10  continue 

numatrmax  *  massburned/(density*denominator) 

do  11,  i  *  0,  (numdivisions-1) ,1 

numberdensity(i)  *  munnormalized(i)  *  numatrmax 

*  /cloudvolume 

11  continue 

open  (1,  file  *  ’population',  status  *  'new') 
do  12,  i  *  0,  (numdivisions-1),  1 
urite(l,*)  numberdensity(i) 

12  continue 
close(l) 

end 

************************ ************* ******41*** ******** ********* 

♦Distribution  finds  the  distribution  of  the  particles. 

♦Called  on  by  part iclepopulat ion .  If  a  new  distribution  is 
♦to  be  studied,  this  is  the  place  to  insert  that 
♦new  distribution. 

♦Inputs:  a  (the  radius  of  the  particles),  step  (increment) 
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*Output:  the  value  of  the  distribution  at  (a)  times  the  step 

*  (area  under  the  curve) 

**************************************************************** 

double  precision  function  distribution(a.step) 

double  precision  a,  step,  distro,  Diam 
Diam  *  2  *  a 

distro  *  (0.051  *  (Diam)**2)  * 

*  exp(  (-1 .89D-02)*((Diam)**2)  ) 

distribution  *  distro  *  step 
end 

**************************************************************** 
★The  following  subroutine  is  the  most  involved  subroutine  in  the 
♦program.  It  may  be  considered  the  heart  of  the  program. 

♦For  a  given  size  of  particles,  it  calculates  the  Riccati-Bessel 
♦functions  needed  to  find  the  Mie  coefficients.  Finds  the 
♦coefficient  of  extinction,  QEXT. 

♦Inputs:  x  *  a/lambda  (particle  radius/wavelength) 

♦  m  *  index  of  refraction 

♦Outputs:  coefficient  of  extention,  QEXT,  mie  coefficients 
♦a  and  b  in  the  file  aandb.dat 

♦Calls  on  subroutine  amplitude,  which  calculates  il  and  i2 
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**************************************************************** 

SUBROUTINE  qAHP(X,M,  QEXT) 

COMPLEX  M,  SY(0:50) ,  SYP(0:50),  Y,  ZETA(0:50) 

COMPLEX  ZETAP(0: 50) ,  A(50),  B(50) 

DOUBLE  PRECISION  LAMBDA,  PI,  X,  SX(0:50),  SXP(0:50) 

DOUBLE  PRECISION  SUMQSCA,  QSCA,  OLDQSCA,  XOFTHETA,  THETA 

DOUBLE  PRECISION  CX(0:50),  CXP(0 :50) ,QEXT,  OLDQEXT 

DOUBLE  PRECISION  SUMQEXT ,  TOLERANCE 

INTEGER  N.MAXITS, LIMIT 

INTRINSIC  SIN,  COS,  ABS 

PARAMETER  (PI  -  3.1415926356) 

PARAMETER  (LIMIT  *  50) 

PARAMETER  (TOLERANCE  »  0.00001) 

EXTERNAL  INITIALRICCATI ,  INITIALRICCATI1 

EXTERNAL  QEXTING,  DERIVATIVES,  AANDB 

EXTERNAL  Q SC ATT, AMPLITUDE,  RICCATIBESSEL 


X  «  2*PI*X 
Y  -  M*X 
SUMQEXT  »  0 
SUMQSCA  -  0 
QEXT  *  0 
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QSCA  «  0 
MAXITS  »  10 

CALL  INITIALRICCATICX,  Y, 

*  CX(0),  SX(0),  SY (0) ,  ZETA(O) , 

*  CXP(O),  SXP(O),  SYP(O) ,  ZETAP(O) ) 

★Comment:  the  above  routine  sets  up  the  initial  Riccati-Bessel 
*  function  and  its  derivative. 

N  »  1 

CALL  INITI ALRICCATI 1 (X , Y ,  CX(l) ,SX(1) ,SY(1) ,ZETA(1)) 

CALL  DERIVATIVES (SX(N-l) ,SX(N) ,CX(N-1) ,CX(N) ,SY(N-1) , 

*  SY (N) ,N,X,Y,  CXP(N) ,SXP(N) ,SYP(N) ,ZETAP(N)) 

1492  DO  10,  N  *  2,  MAXITS 

CALL  RICCATIBESSEL (N , X , Y , SX (N- 1 ) ,SX(N-2) ,CX(N-1) , 

*  CX(N-2) ,SY(N-1) ,SY(N-2) , 

*  CX  (N)  ,  SX  (N)  ,  SY  (N)  ,  ZETA  (N)  ) 


CALL  DERIVATIVES (SX(N-l) ,SX(N) ,CX(N-1) ,CX(N) ,SY(N-1) , 
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*  SY(N),N,X,Y,  CXP(N) ,SXP(N) ,SYP(N) ,ZETAP(N) ) 

10  CONTINUE 


DO  11,  N  *  1,  MAXITS 

CALL  AANDBCSX(N) ,SXP(N) ,SY(N) ,SYP(N) ,ZETA(N) , 
*  ZETAP(N) ,M,  A(N) ,B(N)) 

11  CONTINUE 


DO  12,  N  -  1,  MAXITS 
OLDQEXT  -  QEXT 

CALL  QEXTING(A(N),  B(N) ,  SUMQEXT,  X,  N, 

*  qEXT,  SUMQEXT) 

if  (Cabs (qext  -  oldqext)) .le. (tolerance))  then 

GO  TO  13 
END  IF 

12  CONTINUE 

IF  (MAXITS. Eq. LIMIT)  THEN 
GO  TO  21 
END  IF 


♦COMMENT :  If  MAXITS  has  not  reached  the  limit,  then  increment 
*  by  10  and  go  through  the  process  again. 
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MAXITS  *  MAXITS  +  10 
GO  TO  1492 

21  WRITE (*, 2000) qEXT 

GO  TO  14 

13  WRITE(*, 2010) QEXT 

14  OPEN  (3 .FILE  *  ’ AANDB ’ ,  STATUS  *  'NEW') 

DO  20,  N  *  1,  MAXITS 

WRITE  (3,*)  REAL (A (N)),  AIMAG(A(N) ) , 

*  REAL(B(N)) ,  AIMAG(B(N)) 

20  CONTINUE 

CLOSE  (3) 

CALL  AMPLITUDE (MAXITS) 

2000  FORMAT ( 1 X, 'Q ext  not  converged  after  LIMIT  number  of 

*  terms.  So  far,  qext  *  '.E12.4) 

2010  FORMAT (IX, 'qext  converged  to  '.E12.4) 


END 


*********************************************************** 


*  Subroutines  for  QAMP  follow 

*  Most  of  the  "set  up"  subroutines  axe  self -explanatory . 

*  For  instance,  the  subroutine  initialriccati  finds  the 

*  initial  values  for  the  riccat-bessel  functions,  while  the 

*  subroutine  derivatives  finds  the  derivatives  of  the 

*  riccati-bessel  functions. 

*********************************************************** 

SUBROUTINE  INITIALRICCATI (X,Y, 

*  CX,  SX,  SY,  ZETA, 

*  CXP,  SXP,  SYP,  ZETAP) 


INTRINSIC 

REAL 

COMPLEX 

COMPLEX 


SIN,  COS 

X,  CX,  SX,  CXP,  SXP 

Y,  SY,  ZETA 
SYP,  ZETAP 


SX  *  SIN(X) 

CX  -  COS(X) 

SY  -  SIN(Y) 

ZETA  *  CMPLX(SX.CX) 
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SXP  -  COS(X) 

CXP  *  -SIN(X) 

SYP  *  COS(Y) 

ZETAP  -  CMPLX ( SXP , CXP ) 

END 


SUBROUTINE  INITIALRICCATI 1 (X ,  Y, 

CX,  SX,  SY,  ZETA) 


REAL  X,  CX,  SX 

COMPLEX  Y,  SY,  ZETA 

INTRINSIC  SIN,  COS 

SX  *  (SIN(X) ) /X  -  COS(X) 

CX  -  (COS(X))/X  +  SIN(X) 

SY  *  (SIN(Y))/Y  -  COS(Y) 
ZETA  «  CMPLX (SX,CX) 

END 


* 


SUBROUTINE  RICCATIBESSEL (N , X , Y , SXNM , SXNMM , CXNM , CXNMM , 

SYNM.SYNMM,  CXN.SXN.SYN.ZETAN) 
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REAL  X , SXNM , SXNMM , CXNM , CXNMM , CXN , SXN 
COMPLEX  ZETAN , Y , SYN , SYNM , SYNMM 
INTEGER  N 

SXN  -  ((2*N-1)/X)*SXNM  -  SXNMM 
CXN  -  ((2*N-1)/X)*CXNM  -  CXNMM 
SYN  «  ((2*N-1)/Y)*SYNM  -  SYNMM 
ZETAN  *  CMPLX(SXN.CXN) 

END 


SUBROUTINE  AANDB (SX ,  SXP,  SY,  SYP,  ZETA,  ZETAP,  MM, 
*  A,  B) 

REAL  SX,  SXP 

COMPLEX  SY,  SYP,  ZETA,  ZETAP,  MM,  A,  B 

A  «  (SYP*SX  -  MM*SY*SXP)/(SYP*ZETA  -  MM* SY* ZETAP) 

B  *  (MM*SYP*SX  -  SY*SXP) / (MM*SYP*ZETA  -  SY*ZETAP) 

END 


* 


SUBROUTINE  QEXTING(A,  B,  TEMPSUM,  XTEMP,  NTEMP, 

QTEMP,  TSUM) 


REAL 


TSUM,  TEMPSUM,  QTEMP,  XTEMP 
COMPLEX  A,  B 

INTEGER  NTEMP 

TSUM  *  ( (2*NTEMP  +  l)*REAL(A  +  B))  +  TEMPSUM 
qTEMP  =  (2/XTEMP**2) *TSUM 

END 


SUBROUTINE  qSCATT(A,  B,  TEMPSUM,  X,  N, 

qSCA,  SUM) 

REAL  TEMPSUM,  SUM,  qSCA,  X 

COMPLEX  A ,  B 

INTEGER  N 

SUM  *  (2*N  +  1)  *  ((ABS(A))**2  +  (ABS(B))**2)  +  TEMPSUM 
qSCA  *  (2/X**2)*SUM 

END 


SUBROUTINE  DERIVATIVES (SXNM,  SXN,  CXNM,  CXN,  SYNM.SYN, 
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* 


N,  X,  Y,  CXPN ,  SXPN,  SYPN,  ZETAPN) 


REAL 

COMPLEX 

INTEGER 


SXNM,  SXN,  CXNM,  CXN,  X,  CXPN,  SXPN 
ZETAPN,  SYNM,  SYN,  SYPN.Y 
N 


SXPN  *  SXNM  -  (N/X) *SXN 
CXPN  =  CXNM  -  (N/X) *CXN 
SYPN  =  SYNM  -  (N/Y) *SYN 
ZETAPN  =  ClPLX (SXPN, CXPN) 


END 


*  Comment:  Subroutine  amplitude  finds  the  scattering  functions 

*  for  the  perpendicular  and  plane  polorized  scattered  waves, 

*  il  and  i2,  respectively. 

*  Inputs:  maxits.  Also  uses  the  data  files  aandb.dat  and 

*  viewwidth.dat 

*  Output:  il  and  i2  in  the  data  file  ilandi2.dat 
************** 

SUBROUTINE  AMPLITUDE (MAX ITS) 

complex  a(0:50),b(0:50), sone , stwo ,oldsone,oldstwo 
complex  temp 
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double  precision  pprime(0:50) ,pdoubleprime(0 : 50) 

double  precision  tau(0:50),  p(0:50),  viewwidth 

double  precision  t olerance, reala, imaga, realb , imagb 

double  precision  i 1 (0 : 1,0 : 18) , i2(0 : 1 ,0 : 18) ,PI , theta, x 

double  precision  viewby2 , angledegrees ,  angleradians 

double  precision  sine(0: 1 ,0 : 18) 

integer  n.maxits, count, side 

parameter  (tolerance  =  0.01) 

parameter  (PI  =  3.141592654) 

intrinsic  COS,  CMPLX.abs ,SIN 

open  (1  .file®’ viewwidth' .status® 'old' ) 
read  (1 ,*) viewwidth 
close(l) 

open(l ,f ile®’ aandb ’ .status* 'old') 
do  13,  n  ®  1,  maxits 

read(l,*,end  ®  14)reala, imaga, realb, imagb 
a(n)  *  CMPLX( real a, imaga) 

b(n)  *  CMPLX (realb, imagb) 

13  continue 

14  closed,  status  *  'delete') 

open(2,  file  ®  'ilandi2',  status  ®  'new') 
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viewby2  *  viewvidth/2.0 
count  *  -1 

do  1990,  theta  *  0,  180,  10 
count  *  count  +  1 

do  09,  side  *  0,  1,  1 

angledegrees  *  theta  -  vievby2  +  (viewby2  *  side  *  2) 
angleradians  *  angledegrees  *  PI/180.0 
x  =  C0S(PI  -  angleradians) 

p(0)  *  1 

p(l)  *  x 

do  10,  n  *  2,  maxits 

p(n)  *  2*x*p(n-l)  -  p(n-2)  -  (x*p(n-l)  -  p(n-2))/n 

10  continue 

pprime(0)  *  0 
pprime(l)  *  1 
do  11,  n  =  2,  maxits 

pprime(n)  *  n*p(n-l)  +  x*pprime(n-l) 

11  continue 

pdoubleprime(O)  »  0 
pdoubleprime(l)  *  0 
do  111,  n  *  2,  maxits 
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pdoubleprime(n)  *  (2*n  -  l)*pprime(n-l)  + 

*  pdoubleprime(n-2) 

111  continue 

do  12,  n  *  0,  maxits 

tau(n)  *  x*pprime(n)  -  (l-x**2)*pdoubleprime(n) 

12  continue 

sone  *  cmplx(0 .0 ,0 . 0) 
do  15,  n  *  1,  maxits 
oldsone  *  sone 

temp  *  a(n)*pprime(n)  +  b(n)*tau(n) 
sone  *  (  temp  *  (2*n  +l)/(n*(n+l))  )  +  oldsone 
if  ((abs(sone  -  oldsone)) .le. (tolerance))  then 
go  to  16 
end  if 

15  continue 

print*, ’WARNING:  SI  (and  il)  not  converged’ 

16  il (side, count)  *  (abs (sone) )**2 

stwo  *  cmplx(0. 0,0.0) 
do  18,  n  *  1,  maxits 
oldstvo  *  stwo 

temp  »  b(n)*pprime(n)  +  a(n)*tau(n) 
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stvo  *  (  temp  *  (2*n  +l)/(n*(n+l))  )  +  oldstvo 
if  ((abs(stvo  -  oldstvo)) .le. (tolerance))  then 
go  to  19 
end  if 

18  continue 

print  *, 'WARNING:  S2  (and  i2)  not  converged' 

19  i2(side, count)  *  (abs(stwo))»*2 

sine (side, count)  *  ABS(SIN(angleradians) ) 

09  continue 

vrite(2,*)il(0, count) ,i2(0, count) ,3ine(0, count) , 

*  il(l, count) ,i2(l, count) ,sine(l .count) 

1990  continue 

close(2) 

end 

♦♦♦ *************** *** ***** ****** ******************************** 
♦The  above  subroutine  was  the  last  subroutine  in  QAMP.  By  now,* 
♦the  program  has  calculated,  for  a  given  particle  radius, 

♦the  extinction  coefficient,  gamma,  the  amplitudes  il  and  i2.  ♦ 
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♦The  following  subroutine  uses  il  and  i2  to  find  the  total 
♦intensity  of  the  scattered  light  reaching  the  sensor.  This 
♦total  intensity  is  the  integral  of  the  intensity  over  the 
♦solid  angle  the  viewer  "sees."  The  integration  is  accomplished 
♦numerically  (of  course),  in  this  subroutine. 

♦Input:  k,  numberdensity ,  cloudradius 

♦Output:  the  intensity  of  scattered  light  for  one  particular 
*  direction  in  the  data  file  intensity.dat 

**************************************************************** 

subroutine  intensitysub(k .numberdensity , cloudradius) 

integer  count, side 

double  precision  i 1 (0 : 1,0 : 18) , i2(0 : 1,0 : 18) 
double  precision  sine(0: 1 ,0 : 18) ,  numberdensity 
double  precision  f oftheta(0 : 1) ,  integral,  viewbyi,  PI 
double  precision  intensityCO: 18) ,viewwidth,k 
double  precision  cloudradius ,geomf actor 
parameter  (PI  *  3.1415926356) 

open  (l.file  *  'viewwidth' , status  *  'old') 
read  (1 ,*)viewwidth 
close  (1) 


open  (l,file 


'ilandi2',  status  *  'old') 
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do  10,  count  *  0,  18 

read  (1 ,*)i 1(0, count) , i2(0, count) ,sine(0, count) , 

*  il (1 , count) ,i2(l, count) ,sine(l .count) 

10  continue 

close  (1,  status  *  'delete') 

viewby2  *  viewwidth/2 .0 

*************************************************************** 
♦The  geometric  factor  depends  only  on  the  solid  angle  that 
♦the  detector  "sees."  Thus,  it  is  independent  of  the  sphape 
♦of  the  cloud,  as  long  as  the  detector's  field  of  view  lies 
♦completely  within  the  cloud 

**************************************************************** 


geomfactor  *  (numberdensity*PI*cloudradius)/(k**2) 
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do  11,  count  ■  0,  18 
do  12,  side  *  0,  1 

foftheta(side)*sine(side,count)*(  il (side, count)  + 

i2(side .count)  ) 

continue 


integral  *  (  foftheta(O)  +  foftheta(l)  )  * 

*  (  viewby2  ♦  PI  /  180.0) 
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intensity(count)  *  geomfactor  *  integral 

11  continue 

open  (2, file  *  'intensity',  status  *  'new') 
do  13,  count  *  0,  18 

write  (2 ,*) (intensity(count) ) 

13  continue 

close(2) 
end 
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